rm(list = ls(all = TRUE))
# Estimates the CES parameter #
graphics.off()
library(Rsolnp)
set.seed(20031969)
options(warn=-1)

loss.func.characteristics <-function(params){ 
	beta <- params[1:J]
	sigma <- params[J+1]
	P.index <- apply((beta^sigma)*(Pi^(1-sigma)),2,sum,na.rm=TRUE)   # JPN slide 23 #
	Z.hat <- (beta^sigma)*((Pi/P.index)^(1-sigma))*(X/P.index)		 # JPN slide 25
	Z.hat[which(is.na(Z.hat))]<-0
	e <- as.vector(Z-Z.hat)
	return(t(e)%*%e)
}

loss.func.products <-function(params){ 
	beta <- params[1:K]
	sigma <- params[K+1]
	P.index <- apply((beta^sigma)*(P^(1-sigma)),2,sum,na.rm=TRUE)   
	Q.hat <- (beta^sigma)*((P/P.index)^(1-sigma))*(X/P.index)		
	Q.hat[which(is.na(Q.hat))]<-0
	e <- as.vector(Q-Q.hat)
	return(t(e)%*%e)
}

# ------------------------------------------------- #
# Comment out to use the desired product definition #
# ------------------------------------------------- #

load(file = "results36.RData")         # Broad product definition
load(file = "results171.RData")		# Intermediate product definition
load(file = "results353.RData")		# Narrow product definition


Q[which(is.na(Q))]<-0
T <- dim(P)[2]
J <- dim(Z)[1]
K <- dim(Q)[1]

Pi <- matrix(NA,J,T)
for (t in 1:T){
	Pi[,t] <- apply(Pi.JK[[t]],1,mean,na.rm=TRUE) 
}

colnames(Pi)<-colnames(Z)
rownames(Pi)<-rownames(Z)


# ------------------------------------------------------------------------------ #
# 1. Estimate the CES parameter(s) for the preferences-for-characteristics model #
# ------------------------------------------------------------------------------ #
X<-apply(Pi*Z,2,sum,na.rm=TRUE)
params.0 <- c(apply(Z,1,mean),2)
Upper 	 <- rep(Inf,J+1)
Lower    <- c(rep(0,J),1)
sink("/dev/null")
	CES 	 <- solnp(pars=params.0, loss.func.characteristics, eqfun = NULL, eqB = NULL, ineqfun =  NULL, ineqLB=NULL, ineqUB=NULL, LB=Lower, UB=Upper)
sink()

print("Characteristics")
print(CES$pars[J+1])

# Repeat for each Jack-knife re-samples to get the variance #


	sigma<-NULL
	for (i in 1:K){
		Pi <- matrix(NA,J,T)
		for (t in 1:30){
			Pi[,t] <- Pi.JK[[t]][,i]
		}
		
sink("/dev/null")
		CES <- solnp(pars=params.0, loss.func.characteristics, eqfun = NULL, eqB = NULL, ineqfun =  NULL, ineqLB=NULL, ineqUB=NULL, LB=Lower, UB=Upper)
sink()
		sigma <- c(sigma,CES$pars[J+1])
	}

print(sd(sigma))

print("Goods")

# -------------------------------------------------------------------- #
# 2. Estimate the CES parameter(s) for the preferences-for-goods model #
# -------------------------------------------------------------------- #

X<-apply(P*Q,2,sum,na.rm=TRUE)
params.0 <- c(apply(Q,1,mean),3)
Upper 	 <- rep(Inf,K+1)
Lower    <- c(rep(0,K),1)
sink("/dev/null")
	CES <- solnp(pars=params.0, loss.func.products, eqfun = NULL, eqB = NULL, ineqfun =  NULL, ineqLB=NULL, ineqUB=NULL, LB=Lower, UB=Upper)
sink()

print(CES$pars[K+1])

Price <- P
Quantity <- Q
sigma <- NULL
# Bootstrap the standard errors on sigma #
# sample size is 25%

for (i in 1:200){
	sel <- sample(1:K,round(0.25*K),replace=TRUE)
	P <- Price[sel,]
	Q <- Quantity[sel,]
	X<-apply(P*Q,2,sum,na.rm=TRUE)
sink("/dev/null")		
	CES <- solnp(pars=params.0, loss.func.products, eqfun = NULL, eqB = NULL, ineqfun =  NULL, ineqLB=NULL, ineqUB=NULL, LB=Lower, UB=Upper)
sink()		
	sigma <- c(sigma,CES$pars[K+1])
}

print(sd(sigma))